#ifndef MATHENGINE_H
#define MATHENGINE_H

#include <cmath>
#include <ostream>
#include <cstdio>
#include <type_traits>
#include <algorithm>
#include <iterator>
#include <cassert>
#include <random>

namespace mathengine {

#define GENERAL_MATRIX 1
#define SQUARE_MATRIX  2
#define MATRIX_TYPE    GENERAL_MATRIX

template<typename T>
class MathObject {

public:
    MathObject() {
#if __cplusplus >= 201703L
        if (std::is_same_v<T, float>) {
            useEpsilon = true;
            epsilon = kFloatEpsilon;
        }
        else if (std::is_same_v<T, double>) {
            useEpsilon = true;
            epsilon = kDoubleEpsilon;
        }
        else {
            useEpsilon = false;
            epsilon = 0;
        }
#else

#endif
    }

    T getEpsilon() {
        if (useEpsilon) return epsilon;
        else return 0;
    }

protected:

    static constexpr float kFloatEpsilon = 1e-4;
    static constexpr double kDoubleEpsilon = 1e-8;
    bool useEpsilon;
    T epsilon;
};

template<typename T, size_t N = 2>
class Vec : public MathObject<T> {
public:
    Vec() : dim(N) {
        std::fill(data, data + N, 0);
    }
    explicit Vec(T xx) : dim(N) {
        std::fill(data, data + N, xx);
    }
    Vec(T list[N]) : dim(N) {
        std::copy_n(list, N, data);
    }
    Vec(std::initializer_list<T> list) : dim(N) {
        int copy_size = N > list.size() ? list.size() : N;
        std::copy_n(list.begin(), copy_size, data);
        if (N > copy_size) {
            std::fill(data + copy_size, data + N, 0);
        }
    }
    Vec(const Vec& v) : dim(N) {
        assert(N == v.dim);
        std::copy(v.data, v.data + N, data);
    }
    template<typename U>
    Vec(const Vec<U, N>& v) : dim(N) {
        for (int i = 0; i < N; i++) {
            data[i] = static_cast<T>(v.data[i]);
        }
    }

    Vec& operator=(const Vec& v) {
        assert(N == v.dim);
        if (this != &v) {
            std::copy(v.data, v.data + N, data);
        }
        return *this;
    }
    template<typename U>
    Vec& operator=(const Vec<U, N>& v) {
        assert(N == v.dim);
        for (int i = 0; i < N; i++) {
            data[i] = static_cast<T>(v.data[i]);
        }
        return *this;
    }

    Vec operator+(const Vec& v) const {
        assert(N == v.dim);
        Vec res;
        for (int i = 0; i < N; i++) {
            res.data[i] = data[i] + v.data[i];
        }
        return res;
    }
    Vec operator-(const Vec& v) const {
        assert(N == v.dim);

        Vec res;
        T nums[N];
        for (int i = 0; i < N; i++) {
            res.data[i] = data[i] - v.data[i];
        }
        return res;
    }

    // scalar multiply
    Vec operator*(float s) const {
        Vec res;
        for (int i = 0; i < N; i++) {
            res.data[i] = data[i] * s;
        }
        return res;
    }

    // scalar division, be care s is 0
    Vec operator/(float s) const {
        Vec res;
        for (int i = 0; i < N; i++) {
            res.data[i] = data[i] / s;
        }
        return res;
    }

    Vec& operator+=(const Vec& v) {
        assert(N == v.dim);

        for (int i = 0; i < N; i++) {
            data[i] = data[i] + v.data[i];
        }
        return *this;
    }
    Vec& operator-=(const Vec& v) {
        assert(N == v.dim);

        for (int i = 0; i < N; i++) {
            data[i] = data[i] - v.data[i];
        }
        return *this;
    }
    Vec& operator*=(T s) {
        for (int i = 0; i < N; i++) {
            data[i] = data[i] * s;
        }
        return *this;
    }
    Vec& operator/=(T s) {
        for (int i = 0; i < N; i++) {
            data[i] = data[i] / s;
        }
        return *this;
    }

    Vec operator-() const {
        Vec res;
        for (int i = 0; i < N; i++) {
            res.data[i] = -data[i];
        }
        return res;
    }

    T operator*(const Vec& r) const {
        return dotProduct(r);
    }

    bool operator==(const Vec& v) const {
        if (this->useEpsilon) {
            for (int i = 0; i < N; i++) {
                if (std::abs(data[i] - v.data[i]) > this->epsilon)
                    return false;
            }
        }
        else {
            for (int i = 0; i < N; i++) {
                if (data[i] != v.data[i])
                    return false;
            }
        }

        return true;
    }

    bool operator!=(const Vec& v) const {
        return !operator==(v);
    }

    T dotProduct(const Vec& r) const {
        assert(N == r.dim);
        T sum = 0;
        for (int i = 0; i < N; i++) {
            sum += data[i] * r.data[i];
        }
        return sum;
    }

    Vec<T, 3> crossProduct(const Vec<T, 3>& r) const {
        return Vec<T, 3>({y() * r.z() - z() * r.y(),
                          z() * r.x() - x() * r.z(),
                          x() * r.y() - y() * r.x()});
    }

    T norm() const {
        T sum = 0;
        for (int i = 0; i < N; i++) {
            sum += data[i] * data[i];
        }
        return sum;
    }
    T length() const {
        return std::sqrt(norm());
    }
    Vec& normalize() {
        T n = norm();
        if (n > 0) {
            this->operator/=(std::sqrt(n));
        }
        return *this;
    }

    Vec normalized() const {
        Vec res(*this);
        res.normalize();
        return res;
    }

    T operator[](int i) const {
        assert(i >= 0 && i < N);
        return data[i];
    }
    T& operator[](int i) {
        assert(i >= 0 && i < N);
        return data[i];
    }

    int getDim() const { return dim; }

    T x() const {
        if (N > 0)
            return data[0];
        else
            return std::numeric_limits<T>::infinity();
    }
    T y() const  {
        if (N > 1)
            return data[1];
        else
            return std::numeric_limits<T>::infinity();
    }
    T z() const  {
        if (N > 2)
            return data[2];
        else
            return std::numeric_limits<T>::infinity();
    }
    T w() const {
        if (N > 3)
            return data[3];
        else
            return std::numeric_limits<T>::infinity();
    }

    static Vec zero() {
        return Vec(static_cast<T>(0));
    }

    static Vec random(T min, T max) {
        std::uniform_real_distribution<T> u(min, max);
        std::default_random_engine e;
        Vec res;

        for (int i = 0; i < N; i++) {
            res.data[i] = u(e);
        }

        return res;
    }


    friend std::ostream& operator<<(std::ostream& os, const Vec<T, N>& v) {
        for (int i = 0; i < N - 1; i++) {
            os << v.data[i] << ", ";
        }
        return os << v.data[N-1];
    }

    friend Vec<T, N> operator*(const T& s, const Vec<T, N>& v) {
        return v * s;
    }

public:
    T data[N];
    size_t dim = N;
};


typedef Vec<float, 2>  Vec2f;
typedef Vec<double, 2> Vec2d;
typedef Vec<int, 2>    Vec2i;

typedef Vec<float, 3>  Vec3f;
typedef Vec<double, 3> Vec3d;
typedef Vec<int, 3>    Vec3i;

typedef Vec<float, 4>  Vec4f;
typedef Vec<double, 4> Vec4d;
typedef Vec<int, 4>    Vec4i;


#if (MATRIX_TYPE == GENERAL_MATRIX)

template<typename T, size_t ROWDIM, size_t COLDIM>
class Mat : public MathObject<T> {
public:
    Mat() {
        for (int i = 0; i < ROWDIM; i++) {
            std::fill(data_[i], data_[i] + COLDIM, 0);
        }
    }

    Mat(T xx) {
        for (int i = 0; i < ROWDIM; i++) {
            std::fill(data_[i], data_[i] + COLDIM, xx);
        }
    }

    Mat(T list[ROWDIM * COLDIM]) {
        for (int i = 0; i < ROWDIM; i++) {
            std::copy_n(std::begin(list[0]) + i * COLDIM, COLDIM, std::begin(data_[i]));
        }
    }

    Mat(std::initializer_list<T> list) {
        auto it = list.begin();
        for (int i = 0, count = 0; i < ROWDIM; i++) {
            for (int j = 0; j < COLDIM; j++) {
                data_[i][j] = count < list.size() ? (*it++): 0;
                count++;
            }
        }
    }

    Mat(const Mat& v) {
        assert(v.rowDim_ == rowDim_);
        assert(v.colDim_ == colDim_);

        if (this != &v) {
            for (int i = 0; i < ROWDIM; i++) {
                for (int j = 0; j < COLDIM; j++) {
                    data_[i][j] = v.data_[i][j];
                }
            }
        }

#if 0
        for (int i = 0; i < ROWDIM; i++) {
            std::copy(std::begin(v.data_[i]), std::end(v.data_[i]), std::begin(data_[i]));
        }
#endif
    }

    Mat& operator=(const Mat &v) {
        assert(v.rowDim_ == rowDim_);
        assert(v.colDim_ == colDim_);

        if (this != &v) {
            for (int i = 0; i < ROWDIM; i++) {
                for (int j = 0; j < COLDIM; j++) {
                    data_[i][j] = v.data_[i][j];
                }
            }
        }
#if 0
        for (int i = 0; i < ROWDIM; i++) {
            std::copy(std::begin(v.data_[i]), std::end(v.data_[i]), std::begin(data_[i]));
        }
#endif
        return *this;
    }

    Mat operator+(const Mat& v) const {
        assert(v.rowDim_ == rowDim_);
        assert(v.colDim_ == colDim_);

        Mat res;
        for (int i = 0; i < ROWDIM; i++) {
            for (int j = 0; j < COLDIM; j++) {
                res.data_[i][j] = data_[i][j] + v.data_[i][j];
            }
        }
        return res;
    }

    Mat operator-(const Mat &v) const {
        Mat res;
        for (int i = 0; i < ROWDIM; i++) {
            for (int j = 0; j < COLDIM; j++) {
                res.data_[i][j] = data_[i][j] - v.data_[i][j];
            }
        }
        return res;
    }

    // scalar multiply
    Mat operator*(float s) const {
        Mat res;
        for (int i = 0; i < ROWDIM; i++) {
            for (int j = 0; j < COLDIM; j++) {
                res.data_[i][j] = data_[i][j] * s;
            }
        }
        return res;
    }

    // scalar division, no checking if s == 0
    Mat operator/(float s) const {
        Mat res;
        for (int i = 0; i < ROWDIM; i++) {
            for (int j = 0; j < COLDIM; j++) {
                res.data_[i][j] = data_[i][j] / s;
            }
        }
        return res;
    }

    Vec<T, ROWDIM> operator*(const Vec<T, COLDIM> &v) const {
        Vec<T, ROWDIM> res{0};
        for (int i = 0; i < ROWDIM; i++) {
            for (int j = 0; j < COLDIM; j++) {
                res[i] += data_[i][j] * v[j];
            }
        }
        return res;
    }

    // Mat<T, ROWDIM, COLDIM> * Mat<T, ROWDIM2, COLDIM2>
    // => Mat<T, ROWDIM, COLDIM2>
    template<size_t ROWDIM2, size_t COLDIM2>
    Mat<T, ROWDIM, COLDIM2> operator*(const Mat<T, ROWDIM2, COLDIM2> &v) const {
        static_assert(ROWDIM2 == COLDIM);
        Mat<T, ROWDIM, COLDIM2> res;
        for (int i = 0; i < ROWDIM; i++) {
            for (int j = 0; j < COLDIM2; j++) {
                T sum = 0;
                for (int k = 0; k < ROWDIM2; k++) {
                    sum += data_[i][k] * v.data_[k][j];
                }
                res.data_[i][j] = sum;
            }
        }
        return res;
    }

    Mat& operator+=(const Mat &v) {
        for (int i = 0; i < COLDIM; i++) {
            for (int j = 0; j < ROWDIM; j++) {
                data_[i][j] += v.data_[i][j];
            }
        }
        return *this;
    }

    Mat& operator-=(const Mat &v) {
        for (int i = 0; i < COLDIM; i++) {
            for (int j = 0; j < ROWDIM; j++) {
                data_[i][j] -= v.data_[i][j];
            }
        }
        return *this;
    }

    Mat& operator*=(float s) {
        for (int i = 0; i < ROWDIM; i++) {
            for (int j = 0; j < COLDIM; j++) {
                data_[i][j] *= s;
            }
        }
        return *this;
    }

    // no checking if s == 0
    Mat& operator/=(float s) {
        for (int i = 0; i < ROWDIM; i++) {
            for (int j = 0; j < COLDIM; j++) {
                data_[i][j] /= s;
            }
        }
        return *this;
    }

    bool operator==(const Mat &v) const {
        if (this->useEpsilon) {
            for (int i = 0; i < ROWDIM; i++) {
                for (int j = 0; j < COLDIM; j++) {
                    if (std::abs(data_[i][j] - v.data_[i][j]) > this->epsilon)
                        return false;
                }

            }
        }
        else {
            for (int i = 0; i < ROWDIM; i++) {
                for (int j = 0; j < COLDIM; j++) {
                    if (data_[i][j] != v.data_[i][j])
                        return false;
                }
            }
        }
        return true;
    }

    T getValue(int rowIdx, int colIdx) const {
        assert(rowIdx < ROWDIM && rowIdx >= 0);
        assert(colIdx < COLDIM && colIdx >= 0);

        return data_[rowIdx][colIdx];
    }

    T& getValue(int rowIdx, int colIdx) {
        assert(rowIdx < ROWDIM && rowIdx >= 0);
        assert(colIdx < COLDIM && colIdx >= 0);

        return data_[rowIdx][colIdx];
    }

    Vec<T, COLDIM> getRow(int rowIdx) const {
        assert(rowIdx < ROWDIM && rowIdx >= 0);
        Vec<T, COLDIM> row;
        for (int i = 0; i < COLDIM; i++) {
            row[i] = data_[rowIdx][i];
        }
        return row;
    }
    Vec<T, ROWDIM> getCol(int colIdx) const {
        assert(colIdx < COLDIM && colIdx >= 0);
        Vec<T, ROWDIM> col;
        for (int i = 0; i < ROWDIM; i++) {
            col[i] = data_[i][colIdx];
        }
        return col;
    }

    void setValue(int rowIdx, int colIdx, T v) {
        assert(rowIdx < ROWDIM && rowIdx >= 0);
        assert(colIdx < COLDIM && colIdx >= 0);

        data_[rowIdx][colIdx] = v;
    }

    void setRow(int idx, const Vec<T, COLDIM> &v) {
        assert(idx < ROWDIM && idx >= 0);
        for (int i = 0; i < COLDIM; i++) {
            data_[idx][i] = v[i];
        }
    }

    void setCol(int idx, const Vec<T, ROWDIM> &v) {
        assert(idx < COLDIM && idx >= 0);
        for (int i = 0; i < ROWDIM; i++) {
            data_[i][idx] = v[i];
        }
    }

    // transpose other
    Mat<T, COLDIM, ROWDIM> transposed() const {
        Mat<T, COLDIM, ROWDIM> res;
        for (int i = 0; i < COLDIM; i++) {
            for (int j = 0; j < ROWDIM; j++) {
                res.data_[i][j] = data_[j][i];
            }
        }
        return res;
    }

    // transpose itself
    Mat& transpose() {
        assert(ROWDIM == COLDIM);
        for (int i = 0; i < ROWDIM; i++) {
            for (int j = i + 1; j < COLDIM; j++) {
                std::swap(data_[i][j], data_[j][i]);
            }
        }
        return *this;
    }

    Mat& invert() {
        static_assert(ROWDIM == COLDIM, "the Matrix is not a Square Matrix");
        Mat res = inverse();
        *this = res;
        return *this;
    }

    // A^-1 = A^* / |A|
    Mat inverse() const {
        static_assert(ROWDIM == COLDIM, "the Matrix is not a Square Matrix");
        return (adjoint() * 1.0 / det());
    }

    Mat invert_transpose() const {
        return inverse().transpose();
    }

    // A^* =
    // (A11 A21 ... An1
    //  A12 A22 ... An2
    //  ...
    //  A1n A2n ... Ann)
    //  Aij is the Algebraic Cofactor of the (i,j) matrix A
    Mat adjoint() const {
        Mat res{static_cast<T>(0)};
        for (int i = 0; i < ROWDIM; i++) {
            for (int j = 0; j < COLDIM; j++) {
                res.data_[i][j] = cofactor(j, i);
            }
        }
        return res;
    }

    /// Laplace theorem
    /// |A| = ai1*Ai1 + ai2*Ai2 +...+ ain*Ain
    ///
    /// \see https://www.cnblogs.com/fortunely/p/18076029
    T det() const {
        static_assert(ROWDIM == COLDIM, "det needs Squre Matrix");
        static_assert(ROWDIM * COLDIM <= 1000 * 1000, "Matrix size is too large");

        if (ROWDIM == 1) return data_[0][0];
        if (ROWDIM == 2) return (data_[0][0] * data_[1][1] - data_[0][1] * data_[1][0]);
        T res = 0;
        if (this->rowDim_ > 2) {
            for (int j = 0; j < ROWDIM; j++) {
                //                res += data_[0][j] * cofactor(0, j);
                res += data_[0][j] * ((0 + j ) % 2 == 0 ? 1 : -1) * minorMat(0, j).det();
            }
        }
        return res;
    }

    // Aij = (-1) ^ {i+j} * Mij
    T cofactor(int row, int col) const {
        int sign = ((row + col) % 2 == 0) ? 1 : -1;
        return sign * minorDet(row, col);
    }

    // Mij, ignore A(row; ), A( ;col)
    Mat<T, ROWDIM-1, COLDIM-1> minorMat(int row, int col) const {
        static_assert(ROWDIM == COLDIM, "The matrix must be Squre Matrix");
        Mat<T, ROWDIM-1, COLDIM-1> res;
        for (int i = 0; i < ROWDIM-1; i++) {
            for (int j = 0; j < COLDIM-1; j++) {
                res.data_[i][j] = data_[i < row ? i : i+1][j < col ? j : j+1];
            }
        }
        return res;
    }

    // |Mij|: det of matrix Mij
    T minorDet(int row, int col) const {
        static_assert(ROWDIM == COLDIM, "The matrix must be Squre Matrix");
        static_assert(ROWDIM > 1, "Matrix size must be greater than 1 for minorDet");

        if (ROWDIM == 1) return data_[0][0];
        return minorMat(row, col).det();
    }

    int getRowDim() { return rowDim_; }
    int getColDim() { return colDim_; }

    static Mat zero() {
        return Mat(static_cast<T>(0));
    }
    static Mat identity() {
        Mat res{0};

        for (int i = 0; i < ROWDIM; i++) {
            res.data_[i][i] = static_cast<T>(1);
        }
        return res;
    }

    static Mat random(T min, T max) {
        std::uniform_real_distribution<T> u(min, max);
        std::default_random_engine e;
        Mat res;
        for (int i = 0; i < ROWDIM; i++) {
            for (int j = 0; j < COLDIM; j++) {
                res.data_[i][j] = u(e);
            }
        }
        return res;
    }

    friend std::ostream& operator<<(std::ostream &os, const Mat &v) {
        for (int i = 0; i < ROWDIM; i++) {
            for (int j = 0; j < COLDIM; ++j) {
                os << v.data_[i][j];
                if (j < COLDIM - 1)
                    os << ", ";
            }
            os << "\n";
        }
        return os;
    }

public:
    T data_[ROWDIM][COLDIM];
    size_t rowDim_ = ROWDIM;
    size_t colDim_ = COLDIM;
};

// Partial Specialization, ROWDIM = COLDIM = 1
template<typename T>
class Mat<T, 1, 1> : public MathObject<T> {
public:
    T data_[1][1];

    T det() const {
        return data_[0][0];
    }
};

// Partial Specialization, ROWDIM = COLDIM = 2
template<typename T>
class Mat<T, 2, 2> : public MathObject<T> {
public:
    T data_[2][2];

    T det() const {
        return data_[0][0] * data_[1][1] - data_[0][1] * data_[1][0];
    }

    Mat<T, 1, 1> minorMatrix(int row, int col) const {
        Mat<T, 1, 1> res;
        res.data_[0][0] = data_[row < 1 ? 0 : 1][col < 1 ? 0 : 1];
        return res;
    }

    T minorDet(int row, int col) const {
        return minorMatrix(row, col).det();
    }
};

typedef Mat<float, 2, 2>  Matrix22f;
typedef Mat<double, 2, 2> Matrix22d;
typedef Mat<int, 2, 2>    Matrix22i;

typedef Mat<float, 3, 3>  Matrix33f;
typedef Mat<double, 3, 3> Matrix33d;
typedef Mat<int, 3, 3>    Matrix33i;

typedef Mat<float, 4, 4>  Matrix44f;
typedef Mat<double, 4,4 > Matrix44d;
typedef Mat<int,4, 4>     Matrix44i;

#elif MATRIX_TYPE == SQUARE_MATRIX
template<typename T, size_t N = 2>
class Matrix : public MathObject<T> {
public:

    Matrix() {
        for (int i = 0; i < N; i++) {
            std::fill(m[i], m[i] + N, 0);
        }
    }
    Matrix(T xx) {
        for (int i = 0; i < N; i++) {
            std::fill(m[i], m[i] + N, xx);
        }
    }
    Matrix(T list[N * N]) {
        for (int i = 0; i < N; i++) {
            std::copy_n(std::begin(list[0])+ i * N, N, std::begin(m[i]));
        }
    }
    Matrix(T list[N][N]) {
        static_assert(N > 0 && N <= 1000, "Matrix size is too large or invalid");
        for (int i = 0; i < N; i++) {
            std::copy_n(std::begin(list[i]), N, std::begin(m[i]));
        }
    }
    Matrix(std::initializer_list<T> list) {
        auto it = list.begin();
        for (int i = 0, count = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                if (count < list.size())
                    m[i][j] = *it++;
                else
                    m[i][j] = 0;
                count++;
            }
        }
    }

    Matrix(const Matrix& v) {
        assert(v.dim == N);

        for (int i = 0; i < N; i++) {
            std::copy(v.m[i], v.m[i] + N, m[i]);
        }
    }

    Matrix& operator=(const Matrix& v) {
        assert(v.dim == N);

        for (int i = 0; i < N; i++) {
            std::copy(std::begin(v.m[i]), std::end(v.m[i]), std::begin(m[i]));
        }
        return *this;
    }

    Matrix operator+(const Matrix& v) const {
        Matrix res;
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                res.m[i][j] = m[i][j] + v.m[i][j];
            }
        }
        return res;
    }

    Matrix operator-(const Matrix& v) const {
        Matrix res;
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                res.m[i][j] = m[i][j] - v.m[i][j];
            }
        }
        return res;
    }
    // scalar multiply
    Matrix operator*(T s) const {
        Matrix res;
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                res.m[i][j] = m[i][j] * s;
            }
        }
        return res;
    }

    // scalar division, no checking if s == 0
    Matrix operator/(T s) const {
        Matrix res;
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                res.m[i][j] = m[i][j] / s;
            }
        }
        return res;
    }

    Vec<T, N> operator*(const Vec<T, N>& v) const {
        Vec<T, N> res{0};
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                res[i] += m[i][j] * v[j];
            }
        }
        return res;
    }

    Matrix operator*(const Matrix& v) const {
        Matrix res;
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                T sum = 0;
                for (int k = 0; k < N; k++) {
                    sum += m[i][k] * v.m[k][j];
                }
                res.m[i][j] = sum;
            }
        }
        return res;
    }

    Matrix& operator+=(const Matrix& v) {
        Matrix res;
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                res.m[i][j] = m[i][j] + v.m[i][j];
            }
        }
        *this = res;
        return *this;
    }

    Matrix& operator-=(const Matrix& v) {
        Matrix res;
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                res.m[i][j] = m[i][j] - v.m[i][j];
            }
        }
        *this = res;
        return *this;
    }

    Matrix& operator*=(T s) {
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                m[i][j] = m[i][j] * s;
            }
        }
        return *this;
    }

    Matrix& operator/=(T s) {
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                m[i][j] = m[i][j] / s;
            }
        }
        return *this;
    }

    bool operator==(const Matrix& v) const {
        if (this->useEpsilon) {
            for (int i = 0; i < N; i++) {
                for (int j = 0; j < N; j++) {
                    if (std::abs(m[i][j] - v.m[i][j]) > this->epsilon) {
                        return false;
                    }
                }
            }
        }
        else {
            for (int i = 0; i < N; i++) {
                for (int j = 0; j < N; j++) {
                    if (m[i][j] != v.m[i][j]) {
                        return false;
                    }
                }
            }
        }

        return true;
    }

    void setRow(int idx, const Vec<T, N> &v) {
        assert(idx < N && idx >= 0);
        for (int i = 0; i < N; i++) {
            m[idx][i] = v[i];
        }
    }

    void setColumn(int idx, const Vec<T, N> &v) {
        assert(idx < N && idx >= 0);
        for (int i = 0; i < N; i++) {
            m[i][idx] = v[i];
        }
    }

    /// transpose other
    Matrix transposed() const {
        Matrix res;
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                res.m[i][j] = m[j][i];
            }
        }
        return res;
    }

    /// transpose itself
    Matrix& transpose() {
        for (int i = 0; i < N; i++) {
            for (int j = i + 1; j < N; j++) {
                std::swap(m[i][j], m[j][i]);
            }
        }
        return *this;
    }

    Matrix& invert() {
        Matrix res = inverse();
        *this = res;
        return *this;
    }

    /// A^{-1} = A^* / |A|
    Matrix inverse() const {
        return (adjoint() * 1.0 / det());
    }

    /// A^* =
    /// (A11 A21 ... An1
    ///  A12 A22 ... An2
    ///  ...
    ///  A1n A2n ... Ann)
    ///  Aij is the Algebraic Cofactor of the (i,j) matrix A
    Matrix adjoint() const {
        Matrix res{static_cast<T>(0)};
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                res.m[i][j] = cofactor(j, i);
            }
        }
        return res;
    }


#if 1
    /// Laplace theorem
    /// |A| = ai1*Ai1 + ai2*Ai2 +...+ ain*Ain
    ///
    /// \see https://www.cnblogs.com/fortunely/p/18076029
    T det() const {
        static_assert(N > 0 && N <= 1000, "Matrix size must be greater than 0 for det");

        if (N == 1) return m[0][0];
        if (N == 2) return (m[0][0] * m[1][1] - m[0][1] * m[1][0]);
        T res = 0;
        if (this->dim > 2) {
            for (int j = 0; j < N; j++) {
//                res += m[0][j] * cofactor(0, j);
                res += m[0][j] * ((0 + j ) % 2 == 0 ? 1 : -1) * minorMatrix(0, j).det();
            }
        }
        return res;
    }
#else
    /// LU
    T det() const {
        static_assert(N > 0, "Matrix size must be greater than 0");

        if constexpr (N <= 3) {
            // 对于小矩阵，使用递归展开
            return det();
        } else {
            // 对于大矩阵，使用 LU 分解
            Matrix<T, N> LU = *this;
            T det = 1;

            for (int i = 0; i < N; i++) {
                // 主元不能为 0
                if (LU.m[i][i] == 0) {
                    return 0; // 行列式为 0
                }

                det *= LU.m[i][i]; // 乘以对角线元素

                // 更新下三角和上三角矩阵
                for (int j = i + 1; j < N; j++) {
                    LU.m[j][i] /= LU.m[i][i]; // L 矩阵
                    for (int k = i + 1; k < N; k++) {
                        LU.m[j][k] -= LU.m[j][i] * LU.m[i][k]; // U 矩阵
                    }
                }
            }

            return det;
        }
    }
#endif

    /// Aij = (-1)^{i+j} * Mij
    T cofactor(int row, int col) const {
        int sign = ((row + col) % 2 == 0) ? 1 : -1;
        return sign * minorDet(row, col);
    }

    /// (Mij), ignore A(row; ) and A( ;col)
    Matrix<T, N-1> minorMatrix(int row, int col) const {
        static_assert(N > 1, "Matrix size must be greater than 1 for minorMatrix");

        Matrix<T, N-1> res;
        for (int i = 0; i < N-1; i++) {
            for (int j = 0; j < N-1; j++) {
                res.m[i][j] = m[i < row ? i : i+1][j < col ? j : j+1];
            }
        }
        return res;
    }

    /// |Mij|: determinant of matrix (Mij)
    T minorDet(int row, int col) const {
        static_assert(N > 1, "Matrix size must be greater than 1 for minorDet");

        if (N == 1) return m[0][0];
        return minorMatrix(row, col).det();
    }

    int getDim() { return dim; }


    static Matrix zero() {
        return Matrix(static_cast<T>(0));
    }
    static Matrix identity() {
        Matrix res{0};
        for (int i = 0; i < N; i++) {
            res.m[i][i] = static_cast<T>(1);
        }
        return res;
    }

    static Matrix random(T min, T max) {
        std::uniform_real_distribution<T> u(min, max);
        std::default_random_engine e;
        Matrix res;

        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                res.m[i][j] = u(e);
            }
        }

        return res;
    }

    friend std::ostream& operator<<(std::ostream& os, const Matrix& v) {
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; ++j) {
                os << v.m[i][j];
                if (j < N - 1)
                    os << ", ";
            }
            os << "\n";
        }
        return os;
    }

public:
    T m[N][N];
    int dim = N;
};

// 特化版本：N = 1
template<typename T>
class Matrix<T, 1> : public MathObject<T> {
public:
    T m[1][1];

    T det() const {
        return m[0][0];
    }
};

// 特化版本：N = 2
template<typename T>
class Matrix<T, 2> : public MathObject<T> {
public:
    T m[2][2];

    T det() const {
        return m[0][0] * m[1][1] - m[0][1] * m[1][0];
    }

    Matrix<T, 1> minorMatrix(int row, int col) const {
        Matrix<T, 1> res;
        res.m[0][0] = m[row < 1 ? 0 : 1][col < 1 ? 0 : 1];
        return res;
    }

    T minorDet(int row, int col) const {
        return minorMatrix(row, col).det();
    }
};

typedef Matrix<float, 2>  Matrix22f;
typedef Matrix<double, 2> Matrix22d;
typedef Matrix<int, 2>    Matrix22i;

typedef Matrix<float, 3>  Matrix33f;
typedef Matrix<double, 3> Matrix33d;
typedef Matrix<int, 3>    Matrix33i;

typedef Matrix<float, 4>  Matrix44f;
typedef Matrix<double, 4> Matrix44d;
typedef Matrix<int,4 >    Matrix44i;
#endif

template<typename T, size_t N>
Vec<T, N> operator*(float t, const Vec<T, N> &v) {
    return v * t;
}

template<typename T, size_t N>
Vec<T, N> operator*(const Vec<T, N> &l, const Vec<T, N> &r) {
    return l.dotProduct(r);
}

template<typename T, size_t N>
Vec<T, 3> operator^(const Vec<T, 3> &l, const Vec<T, 3> &r) {
    return l.crossProduct(r);
}

template<typename T, size_t N>
bool operator==(const Vec<T, N>& l, const Vec<T, N>& r) {
    return l.operator==(r);
}

template<typename T, size_t N>
T dotProduct(const Vec<T, N>& l, const Vec<T, N>& r) {
    return l.dotProduct(r);
}

// Vector3(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
template<typename T>
Vec<T, 3> crossProduct(const Vec<T, 3>& l, const Vec<T, 3>& r) {
    return Vec<T, 3>{l.y() * r.z() - l.z() * r.y(),
                     l.z() * r.x() - l.x() * r.z(),
                     l.x() * r.y() - l.y() * r.x()};
    // return l.crossProduct(r);
}

#if 1
template<typename T, size_t N>
inline Vec<T, N> Lerp(const Vec<T, N>& p1, const Vec<T, N>& p2, float t) {
    Vec<T, N> res  = p1 * (1.f - t) + p2 * t;
    return res;
}
#endif

#if MATRIX_TYPE == GENERAL_MATRIX
template<typename T, size_t ROWDIM, size_t COLDIM>
bool operator==(const Mat<T, ROWDIM, COLDIM> &l, const Mat<T, ROWDIM, COLDIM> &r) {
    return l.operator==(r);
}
template<typename T, size_t ROWDIM, size_t COLDIM>
Mat<T, ROWDIM, COLDIM> operator*(T s, const Mat<T, ROWDIM, COLDIM>& mat) {
    return mat.operator*(s);
}

template<typename T, size_t ROWDIM, size_t COLDIM, size_t COLDIM2>
Mat<T, ROWDIM, COLDIM2> multiply(const Mat<T, ROWDIM, COLDIM> &l, const Mat<T, COLDIM, COLDIM> &r) {
    return l.operator*(r);
}

template<typename T, size_t ROWDIM, size_t COLDIM, size_t ROWDIM1, size_t COLDIM2>
Mat<T, ROWDIM, COLDIM2> dotProduct(const Mat<T, ROWDIM, COLDIM> &l, const Mat<T, ROWDIM1, COLDIM2> &r) {
    return l.operator*(r);
}

// Left multiply: row operations, mat * v
template<typename T, size_t ROWDIM, size_t COLDIM>
Vec<T, ROWDIM> operator*(const Mat<T, ROWDIM, COLDIM>& mat, const Vec<T, COLDIM>& v) {
    return mat.operator*(v);
}

#elif MATRIX_TYPE == SQUARE_MATRIX
template<typename T, size_t N>
bool operator==(const Matrix<T, N>& l, const Matrix<T, N>& r) {
    return l.operator==(r);
}

template<typename T, size_t N>
Matrix<T, N> operator*(T s, const Matrix<T, N>& mat) {
    return mat.operator*(s);
}

template<typename T, size_t N>
Matrix<T, N> dotProduct(const Matrix<T, N>& l, const Matrix<T, N>& r) {
    return l.operator*(r);
}

// Left multiply: elementary row operations
template<typename T, size_t N>
Vec<T, N> operator*(const Matrix<T, N>& mat, const Vec<T, N>& v) {
    return mat.operator*(v);
}

template<typename T, size_t N>
Vec<T, N> leftMultiplyVec(const Matrix<T, N>& mat, const Vec<T, N>& v) {
    return mat.operator*(v);
}
#endif

template<typename T, int FROM, int TO>
Vec<T, TO> embed(const Vec<T, FROM>& v, T fill = 1) {
    static_assert(FROM <= TO, "N must be less than or equal to M");
    Vec<T, TO> res;
    for (int i = 0; i < TO; i++) {
        res[i] = i < FROM ? v[i] : fill;
    }
    return res;
}

template<typename T, int FROM, int TO>
Vec<T, TO> proj(const Vec<T, FROM>& v) {
    static_assert(FROM >= TO, "N must be greater than or equal to M");
    Vec<T, TO> res;
    for (int i = 0; i < TO; i++) {
        res[i] = v[i];
    }
    return res;
}

// e.g (v.x, v.y, v.z) -> (v.x, v.y, v.z, 1)
template<typename T, size_t N>
Vec<T, N+1> ToHomogeneousCoord(const Vec<T, N>& v) {
    Vec<T, N+1> res;
    for (int i = 0; i < N; i++) {
        res[i] = v[i];
    }
    res[N] = 1;
    return res;
}

// e.g (v.x, v.y, v.z, v.w) -> (v.x/v.w, v.y/v.w, v.z/v.w)
// require N > 0, v.w != 0
template<typename T, size_t N>
Vec<T, N-1> ToCartesianCoord(const Vec<T, N>& v) {
    assert(N > 0 && v[N-1] != 0);
    Vec<T, N-1> res;
    for (int i = 0; i < N-1; i++) {
        res[i] = v[i] / v[N-1];
    }
    return res;
}

} // namespace mathengine

#endif // MATHENGINE_H
